Proceedings of IDETC/CIE 2007 
ASME 2007 International Design Engineering Technical Conferences & 
Computers and Information in Engineering Conference 
September 4-7, 2007, Las Vegas, USA 


DRAFT DETC2007-3551 1 


INVESTIGATION OF BOUNDARY CONDITIONS FOR FLEXIBLE MULTIBODY 
SPACECRAFT DYNAMICS (DRAFT VERSION) 


John R. MacLean Leslie J. Quiocho* 

An Huynh 

NASA Johnson Space Center 

METECS Automation, Robotics, & Simulation Division 

Houston, Texas 77258 Houston, Texas 77058 

Email: jmaclean@metecs.com, ahuynh@metecs.com Email: leslie.j.quiocho@nasa.gov 


ABSTRACT 

In support of both the Space Shuttle and International Space 
Station programs , a set of generic multibody dynamics algo- 
rithms integrated within the Trick simulation environment have 
addressed the variety of on-orbit manipulator simulation re- 
quirements for engineering analysis, procedures development 
and crew familiarization/training at the NASA Johnson Space 
Center (JSC). Enhancements to these dynamics algorithms are 
now being driven by a new set of Constellation program require- 
ments for flexible multibody spacecraft simulation. One partic- 
ular issue that has been discussed within the NASA community 
is the assumption of cantilever-type flexible body boundary con- 
ditions. This assumption has been commonly utilized within ma- 
nipulator multibody dynamics formulations as it simplifies the 
computation of relative motion for articulated flexible topolo- 
gies. Moreover, its use for modeling of space-based manipula- 
tors such as the Shuttle Remote Manipulator System (SRMS) and 
Space Station Remote Manipulator System (SSRMS) has been ex- 
tensively validated against flight data. For more general flexible 
spacecraft applications, however, the assumption of cantilever- 
type boundary conditions may not be sufficient. This paper de- 
scribes the boundary condition assumptions that were used in 
the original formulation, demonstrates that this formulation can 
be augmented to accommodate systems in which the assump- 
tion of cantilever boundary conditions no longer applies, and 
verifies the approach through comparison with an independent 


model previously validated against experimental hardware test 
data from a spacecraft flexible dynamics emulator. 

NOMENCLATURE 

Tl n Strain energy of body n 
qj The j th generalized DOF of the system 
qj First time derivative of qj 
Qj The j th generalized force 

v * Absolute velocity of point mass dm inside body n 

a* Absolute acceleration of dm inside body n 

q n Flexible DOF of body n 

q n First time derivative of q n 

q n Second time derivative of q n 

M rpn Rigid-rigid mass matrix of body n 

M re ^ n Rigid-elastic mass matrix of body n 

M epn Elastic-rigid mass matrix of body n 

M ee n Elastic-elastic mass matrix of body n 

K ee n Elastic stiffness matrix of body n 

F np Force and moment of body p acting on body n 

F ns Force and moment of body s acting on body n 

r s Spatial vector from joint n to joint s 

E pn External rigid forces and moments acting on body n 

E e:n External elastic forces and moments acting on body n 

B pn Nonlinear rigid inertia force of body n 

B e n Nonlinear elastic inertia force of body n 

S s Mode shape and mode slope of joint s due to body n flex 


Address all correspondence to this author. 


1 


Copyright © 2007 by ASME 



F n Force and moment from body p (previous body to n) acting 
at joint n 

F s Force and moment from body n (previous body to s) acting 
at joint s 

BACKGROUND 

In support of both the Space Shuttle and International Space 
Station programs, a set of generic multibody dynamics algo- 
rithms integrated within the Trick simulation environment have 
addressed the variety of on-orbit manipulator simulation require- 
ments for engineering analysis, procedures development and 
crew familiarization/training at the NASA Johnson Space Cen- 
ter (JSC) [1,2]. Enhancements to these dynamics algorithms 
are now being driven by a new set of Constellation program 
[3] requirements for flexible multibody spacecraft simulation. 
One particular issue that has been discussed within the NASA 
community is the assumption of 4 ‘cantilever- type” flexible body 
boundary conditions [4]. This assumption has been commonly 
utilized within manipulator multibody dynamics formulations as 
it simplifies the computation of relative motion for articulated 
flexible topologies. Moreover, its use for modeling of space- 
based manipulators such as the Shuttle Remote Manipulator Sys- 
tem (SRMS) and Space Station Remote Manipulator System 
(SSRMS) has been extensively validated against flight data. For 
more general flexible spacecraft applications, however, it has 
been pointed out that the assumption of cantilever- type bound- 
ary conditions may not be sufficient [5]. 

The purpose of this paper is to investigate the application 
of the Trick-based multibody dynamics algorithms for a specific 
case in which the assumption of cantilever boundary conditions 
no longer applies. It is shown that the formulation can be re- 
cast to accommodate more general boundary conditions without 
significant modifications or reduction in efficiency. To verify 
the approach, a case study is chosen that includes pinned end 
mass boundary conditions. The reformulated multibody dynam- 
ics model is then correlated against a transfer function model 
which has been previously validated against experimental hard- 
ware data obtained from a spacecraft flexible dynamics emulator. 


MULTIBODY DYNAMICS FORMULATION 

In [1], Modified Lagrange’s Equations were developed to 
support the simulation of space-based manipulators such as the 
SRMS and SSRMS: 


— ; a„dm-\ = Q: 

4j 4j 


( 1 ) 



Figure 1 . BODY N SCHEMATIC 


Mfr^nAn "T" M ten cj fl — F n v S F S -|- F r n -|- B , n (2) 

s^n 


M e r,n,An T~ M een q n -j- K ee n q n — H - E e , n T B en (3) 

s^n 


Equations 2 and 3 represent the core body level elements 
of the Trick-based multibody dynamics package, also known as 
MBDYN. 

Figure 1 depicts a schematic representation of body n. Here 
the inboard frame is given as Frame n while the outboard frame 
is shown as Frame s. F n represents the agregate force and mo- 
ment vector of the previous body acting on body n at Frame n\ 
likewise F s is the force and moment vector of current body acting 
on the next body s at Frame s. Referring back to Equation 3, it is 
important to notice that the S„F n term corresponding to the mode 
shape and mode slope of Frame n is not included. The absence 
of this particular term results from the cantilever-type boundary 
condition discussed in the previous section. 

Extending these equations to the more general boundary 
condition (i.e., eliminating the cantilever- type assumption) re- 
sults in 


Using Equation(l), it was shown that the EOMs in matrix 
form for body n , a single body from within a tree topology, are 


M rr ^ n A n + M re n q n — F n r s^s T E r ^ n -f- B rn 

s^n 


( 4 ) 
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M;,X + M* eeA q n + K een q n = -^S s * T F s + E e , n + B* en (5) 

s^n 

where 


M e r,n — M e r,n ^rr,n^n ( 6 ) 


M e e,n — M ee ,n ^er,n^n $n M re,n H - $n ^rr,n^n (7) 


Ss* T =Ss T -S n T r s (8) 

It is important to note that the original form of the equations 
remains intact, although the definition of matrix expressions has 
now been ‘enhanced’ to eliminate the cantilever-type boundary 
condition limitation. 


CASE STUDY 

In order to test the proposed approach for handling more 
general boundary conditions, a satellite arm flexible dynamics 
emulator was modeled in MBDYN. Whereas the space robotics 
systems previously modeled using MBDYN are characterized by 
stiff links and high ratio gearboxes at the articulations, this emu- 
lator has a very flexible link driven by a brushless direct drive 
motor and hence is a good test case for implementation of a 
pinned boundary condition. In addition, the dynamics system 
is simple enough to be modeled analytically. This section de- 
scribes the emulator, its relevant dynamics characteristics, and 
implementation of the model using MBDYN. 

Experimental Hardware 

As reported in [6], an experimental apparatus was designed 
for the Communications Research Center in Ottawa, Ontario 
Canada to serve as a satellite arm flexible dynamics emula- 
tor. This research was performed in support of the joint United 
States/Canada Mobile Satellite (MS AT) system program (as de- 
picted in by the illustration in Figure 2). The emulator consisted 
of a hinged- free flexible stainless steel beam with a direct-drive 
DC torque motor at the hinged end and a tip mass and position 
sensor at the free end. The beam was free to rotate while the 
tip mass floated on an air bearing table made of granite. The tip 



Figure 2. MS AT SATELLITE 



Figure 3. FLEXIBLE BEAM APPARATUS 


mass was connected to the flexible beam via a pin type bearing. 
The output of the position sensor represented the transverse dis- 
placement of the tip from a reference position in the middle of 
the air bearing table. An additional sensor measuring the angular 
rotation of the motor shaft was also incorporated for experimen- 
tation. The apparatus is shown in Figure 3 while a schematic 
diagram representing the system as seen from above is provided 
in Figure 4. 

Because of the high flexibility of the beam, the absence of 
a gearbox, and the analytically simple dynamics of the system, 
the emulator provides an excellent test case for MBDYN model 
implementation of a pinned boundary condition. 

Unfortunately, as university lab space is frequently in short 
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Sensor 



X (x) —C\ sinh (^j x ^j + C 2 COSI 1 +C 3 sin (^j x ^j + C 4 COS • 

W ; 7 (12) 
The four boundary conditions for the pinned end mass beam 
corresponding to the satellite arm emulator are 


X (0) = 0 


(13) 


Figure 4. FLEXIBLE BEAM SCHEMATIC 

supply, this apparatus has not existed for many years. An analytic 
model, however, was derived in [ 6 ] and was validated against 
the hardware test data, some of which still exists in the form 
of Matlab files. The analytic model, which is described below, 
makes no use of the cantilever-type assumption that restricted the 
MBDYN formulation and hence forms a suitable basis for testing 
the proposed MBDYN modeling approach. 


d 2 X (0) 
dx 2 


= 0 


(14) 


d 2 X (/) 
dx 2 


= 0 


(15) 


Analytic Emulator Model 

This section describes the analytic model that was developed 
in [ 6 ] and validated against the response of the emulator appara- 
tus. As shown below the model makes no cantilever assumption. 
Consider the classical Bemoulli-Euler beam equation [7] 


d 4 y p (x,t) d 2 y p (x,t) 

EI ^^ + 9 ^t^ = 


(9) 


,3 

dx 3 


+ E//X (/) = 0 


(16) 


where // = and E represents the ratio between the tip 
mass and beam mass. For a stable solution we reject the case 
where a < 0 and we consider c = 0 and a > 0 . 

For the first case of a = 0, 


Using separation of variables, Equation (9) can be broken 
into two ordinary differential equations: 


and 


X 0 (x) = A 0 x 


(17) 


El d 4 X (x) 
pX (x) dx 4 


aX (x) = 0 


( 10 ) 


® o (t)=B 0 t+C 0 (18) 


where Ao is an arbitrary constant and Bo, Co are determined 
from the initial conditions. 

For the second case of a > 0, it can be shown that there are 

d 2< & (t) ^ x x an infinite number of solutions O;, i = 1 , 2 , 3. . .00 such that 

J +o<S>(f) = 0 (11) ’ ’ ’ 


where a is a constant. The homogeneous solution for 
Equation(9) is given as 



(19) 
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and 


and 


cot jut = coth jjj + 2 Zfii 


( 20 ) 


yjl ( 1 + ^ + 2 ©/i/ + coth^/) ) 


(27) 


The solutions to Equations (10) and (11) are given as 


Xi(x)=Ai 


sinh (jx) 
sinh jut 


sin (f 4 

sin ^ 


(21) 


and 


Using virtual work techniques outlined in [6], the continuous 
time transfer functions from motor torque to tip position is given 
by 


Ptip (s) 


I 


dXj( 0) 
dx 

P 


Mi) 

s 2 + 2^,-co is + CO / 2 


(28) 


d>/ (t) = Bisin + Q) 


( 22 ) 


and the transfer function from motor torque to hub position 
is given by 


where 


CO/ = 



(23) 


Again, A/ are arbitrary constants and 5/,C/ depend on the 
initial conditions. The values CO/ represent the system modal fre- 
quencies while the functions 2Q (v), or ‘shape’ functions, are the 
system mode shapes. 

It can be shown that these mode shapes, 2Q (x), are orthogo- 
nal with respect to the scalar product 


{Xi(x) | Xj(x)) = j‘r(x)Xj(x)Xj(x)dx (24) 

where 

( l r{x) = [1 H-S/5(jc — /)] (25) 

Jo 

and 8() is the Dirac delta function. 

It is convenient to chose A/ to normalize Xi(x) (i.e., 
(Xi (x) | Xj (x)) = 1). This results in the following values 


Ao 


/ 3 (1 + 3S) 


(26) 


Phub { s ) 


I 


d 2 Xj( 0) 

dx 2 

p(i 2 + 2^C0^ + C0; 2 ) 


(29) 


The contribution of the direct-drive torque motor with am- 
plifier gain K a , torque constant K u winding resistance R and back 
EMF constant Kp is given by 


K a K t K s P tip (s ) 

R + KfKpsPhub (s) 


(30) 


Physical parameters of the system were determined by ex- 
perimentation. They are listed in Table 1 . Recall that E repre- 
sents the ratio of the tip mass to mass of the flexible beam. 

To validate the analytic model implementation, a sweep of 
sinusoidal forcing functions were applied to the system to iden- 
tify the modal frequencies and these frequencies were compared 
to those derived theoretically from the model. This comparison, 
shown in Table 2, demonstrated a match of within 3.4%. 

The model was further verified by measuring the gain and 
the phase delay from input voltage to sensed tip position at a 
range of frequencies. The results shown in Figures 5 and 6 for 
frequencies from 0.1 Hz to 4 Hz demonstrate good agreement 
between analytic and hardware responses. 

One can conclude, therefore, that the analytic emulator 
model, in the form of continuous time transfer functions, has 
been shown to be a good representation of the physical hard- 
ware. Based on the success of this comparison to hardware data, 
the parameters of the analytic model are used in the following 
section to derive input data for the more generalized MBDYN 
formulation which is then compared to the response of the ana- 
lytic model. 
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Table 1 . BEAM CHARACTERISTICS 


1000 


Parameter 

Description 

Value (SI units ) 

E 

Modulus of elasticity 

N/A 

I 

Area moment of inertia 

N/A 

P 

Linear density 

4.09 x 10~ 4 Ib-sec2/in2 



(2.82 kg/m) 

[ei 

V p 

Combined beam properties 

9.3 x 10 3 in2/sec 



(5.974 m2/s) 

l 

Length 

132.75 in 



(3.372 m) 


Mass Ratio 

5.18 



0.1 1 10 100 
rad/s 


Figure 5. COMPARISON OF MAGNITUDE 


Table 2. MEASURED VS. CALCULATED MODAL FREQUENCIES 


i 

co ; (measured) 

jui (E = 5.18) 

co i (calculated) 

1 

5.25 rad/s 

3.171011466 

5.28 rad/s 

2 

21.3 rad/s 

6.298221427 

20.9 rad/s 

3 

48.3 rad/s 

9.434866134 

46.7 rad/s 

4 

86.1 rad/s 

15.71404495 

83.4 rad/s 


MBDYN Model 

To simulate the flexible spacecraft emulator dynamics in 
MBDYN, the matrix quantities M rr ^ n , M re ^ n , M ee ^ n , K ee ^ n , and Sj 
must be derived from the system properties for input to the body 
level equations of motion. 

The beam rigid mass properties, in terms of the M rr ^\ matrix 
is 



Figure 6. COMPARISON OF PHASE 


Evaluating the terms of M re , n for the problem at hand results 
in 


M rr , 1 = 


m\ [ 1 ] —m\c\ 
m\C\ J\ 


(31) 



/ cosh^ — 1 
V sinh w 


cos Hi — 1 
sin jui 


(33) 


The rigid elastic coupling matrix is given as 


and 


M re ,n 


Jo Xi(x)dm 
Jq r x Xi (x) dm 


(32) 


/' 

Jo 


r x X, ( x)dm 


Aj 

sinhp, 


( — cosh m 
\Vi 


y sinh p, 

f4 
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Ai ( l 2 


- cos jm + -J sin /Hi 
sm jut \ m juf 


The elastic mass matrix can be represented by 

M ee ,n = Xi (x)Xj (. X ) dx 

Substituting 2Q from Equation (21) results in 

( ifi \ 

— j 4w IA{A j 

for i 7 ^ j and 


(34) 


(35) 


(36) 
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Figure 7. INPUT TORQUE 

The resulting mode shape and slope matrix is therefore 


/m\ 1 - 35 + 2 Zft (Z/jj + cothjU,) 
V Z / 1 + w + 2 £j/Ui (i/i; + COthpi) 


(37) 


S s (x) = 


X(x) 
a (x) 


(42) 


for i = j 

The stiffness matrix can be found from 


f l d 2 Xi d 2 Xj 

K een = El , ' { dx (38) 

Jo dx z dx z 

resulting in 


K ee , n = El (y ) 4 (39) 

The lateral displacement as a function of distance along the 
beam, x is given as 


The boundary condition at v = 0 corresponds to a nonzero 
S n while the boundary condition at the tip corresponds to S s . To 
implement the model in MBDYN, Equations ( 6 ), (7) and ( 8 ) are 
then used to calculate the required inputs. 

MODEL COMPARISON 

The next step of the process was to correlate the MB- 
DYN implementation of the system model discussed previously 
against the original analytic dynamics model validated against 
the experimental hardware. To perform this model validation, a 
forcing function (see Figure 7) designed to excite the first two 
frequencies was employed and the resulting states for the hub 
and tip position were compared. Comparison plots are provided 
in Figures 8 and 10, while the difference plots are given in Fig- 
ures 9 and 11. Note that an excellent comparison is achieved. 


y(x,t) — X(x)§(t) 

The mode slopes can then be calculated as 

dXi Uj ( cosh (jx) cos(fx) 
a t = —2 = A’J . y + u y 
dx l \ smh/Ji sm jm 


(40) CONCLUDING REMARKS 

To summarize, although originally formulated with the as- 
sumption of cantilever- type modes, it has been demonstrated 
through transformation of body level inputs that MBDYN is ca- 
pable of handling more general boundary conditions using its ex- 
isting mathematical formulation and implementation. Moreover, 

( 41 ) this capability has been verified through a case study by com- 
paring to an analytical model (validated against test hardware) 
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Figure 8. HUB POSITION COMPARISON 


Figure 10. TIP POSITION COMPARISON 



0123456789 10 
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Figure 9. HUB POSITION DIFFERENCE 


Figure 11. TIP POSITION DIFFERENCE 


where the cantilever assumption no longer applies. Finally, it 
is anticipated that the approach described here within this paper 
can be extended to other scenarios involving flexible multibody 
spacecraft simulation with other types of boundary conditions. 
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